%% 现代优化算法Modern optimization algorithm
%% 1.模拟退火Simulated annealing
% 模拟退火算法得益于材料统计力学的研究成果。在高温状态下，粒子能量较高，可以
% 自由移动和重新排列，在低温状态下，粒子能量较低。如果从高温状态下，非常缓慢
% 的降温，粒子就可以在每个温度下达到热平衡。当系统完全被冷却时，最终形成处于
% 低能状态的结晶。
% 下方通过程序说明模拟退火进行过程，程序可见文件：Simulated_annealing.m
clear,clc
% 参数初始化
f=@(x)1/12*sum(x.^4)-1/2*sum(x.^2)+sum(abs(x-4));% 目标函数
lb=-1e2*ones(1,6);ub=1e2*ones(1,6);%变量范围
n=length(lb);% 变量个数
% 模拟退火计算过程
mn=10;% 进行10次，每次初始点不同
mn1=30;% 每个初始点进行30次模拟退火过程
mn2=2e3;% 模拟退火迭代最大次数
n1=4;% 每次模拟退火过程选择随机方向个数
Tm=100;% 最大温度
gama=0.99;% 退火率
result=[];% 存储mn次初始点对应最优结果集
for i=1:mn
    tic
    result1=zeros(mn1,n+1);% 存储mn1次模拟退火过程计算结果
    x0=lb+i/mn*(ub-lb);% 初始点选取
    for i1=1:mn1
        T=Tm;% 温度初始化
        x1=x0;x2=x0;% x1可接受劣解，x2不接受劣解，但迭代时以x1为准
        r2=f(x2);% 计算x2处值
        r1=f(x1);% 计算x1处值
        for i2=1:mn2
            z=-1+2*rand(n1,n);% 随机方向
            for i3=1:n1
                z1=z(i3,:);
                z1=z1/norm(z1);% 方向规范化
                xk=x1+T/Tm*z1.*(ub-lb);
                for i4=1:n
                    if(xk(i4)>ub(i4))
                        xk(i4)=x1(i4)+1/2*(ub(i4)-x1(i4));
                    elseif(xk(i4)<lb(i4))
                        xk(i4)=x1(i4)+1/2*(lb(i4)-x1(i4));
                    end
                end
                r0=f(xk);% 计算xk处值
                if(r0<=r1)
                    x1=xk;
                    r1=r0;
                    if(r0<r2)
                        x2=xk;
                        r2=r0;
                    end
                else
                    p0=exp(-1/T);% 该部分针对不同的目标函数处理方式不
                    % 同，以便x可以接受劣解,如：
                    % p0=exp(-(r0-r1)/T);
                    z2=rand;
                    if(z2>p0)
                        x1=xk;
                        r1=r0;
                    end
                end
            end
            T=gama*T;% 退火
        end
        result1(i1,:)=[r2,x2];
    end
    result=[result;result1];
    toc
end
% 可见，模拟退火算法一般可以找到目标函数的较优解，因此针对具体的复杂优化问题
% 可以首先利用模拟退火算法寻找较优解，再利用梯度线搜索算法或穷举法在较优解附
% 近进行精确搜索。
% 下面本节利用模拟退火寻找TSP问题的较优解，代码及数据来源于CSDN博主
% 「Handsome_Zpp」的原创文章，原文链接为：
% https://blog.csdn.net/qq_44865735/article/details/124100778
% 程序可见文件：SA_TSP.m
% 首先设计计算旅行商距离的函数：
% 计算距离的函数，见文件distance1.m
% 从输入City中从上到下依次距离累加，最后计算最后一个点到起点距离，形成一个圈
% function len=distance1(City,n)
% len=0;
% for i=1:n-1
%     len=len+sqrt(sum((City(i,:)-City(i+1,:)).^2));
% end
% len=len+sqrt(sum((City(n,:)-City(1,:)).^2));
% end
clear,clc                         
% 初始化
load('city.mat');% 导入城市数据
C=city;% 城市数据两列分别为各个城市的x，y轴坐标
n=size(city,1);% 得到TSP问题的规模,即城市数目
T=100*n;% 设置初始温度
L=100;% 马可夫链长度
K=0.99;% 衰减参数（退火速率）
l=1;% 初始化统计迭代次数为1
len(l)=distance1(city,n);% 每次迭代后的路线长度初始化（目标函数值）
% 迭代求解 
figure(1); 
e=1e-2;
len3=distance1(city,n);% 迭代过程最优点，不接受劣解，但不影响迭代
while T>e % 停止迭代温度
    len1=len(l);% 计算旧点值
    for i=1:L % 每一个温度下进行L次方向选取
        % 产生新解       
        p1=floor(1+n*rand());% 随机选择两个不同的城市的坐标
        p2=floor(1+n*rand());
        while p1==p2 % 确保选择的两个城市不一样
            p1=floor(1+n*rand());
            p2=floor(1+n*rand());
        end
        tmp_city=city;% 置换两个不同的城市的坐标
        tmp=tmp_city(p1,:);
        tmp_city(p1,:)=tmp_city(p2,:);
        tmp_city(p2,:)=tmp;
        len2=distance1(tmp_city,n);% 计算新路线的距离 
        delta_e=len2-len1;% 新老距离的差值，相当于能量的差
        % 新路线好于旧路线，用新路线代替旧路线
        if delta_e<0        
            city=tmp_city;
            len1=len2;
            if(len2<len3)
                city1=tmp_city;
                len3=len2;
            end
        else
        % 以概率选择是否接受新解
            if exp(-1/T)>rand()% 该处需结合不同的目标函数进行修改，如：
                % exp(-delta_e/T)>rand()
                city=tmp_city;
                len1=len2;
            end
        end
    end
    l=l+1;% 更新迭代次数
    len(l)=len1;% 计算新路线距离 
    T=T*K;% 更新温度
    % 绘制迭代图
    for i=1:n-1
        plot([city(i,1),city(i+1,1)],[city(i,2),city(i+1,2)],'bo-');
        hold on;
    end
    plot([city(n,1),city(1,1)],[city(n,2),city(1,2)],'ro-');
    title(['最短距离:',num2str(len(l))]);
    hold off;
    pause(0.002);
end
figure(2);
for i=1:n-1
    plot([city(i,1),city(i+1,1)],[city(i,2),city(i+1,2)],'bo-');
    hold on;
end
plot([city(n,1),city(1,1)],[city(n,2),city(1,2)],'ro-');
for i=1:n
    text(C(i,1)+0.5,C(i,2),num2str(i));       %标号  
end
text(city(1,1),city(1,2),'    起点');
text(city(n,1),city(n,2),'    终点');
title('最终路线图');
figure(3);
plot(len)
xlabel('迭代次数')
ylabel('目标函数值')
title('适应度进化曲线')
%% 2.遗传算法Genetic algorithm
% 遗传算法是一种群体搜索技术，根据适者生存原则逐步进化寻找较优解。遗传算法一
% 般包含如下步骤：产生初始群体、求每一个个体适应度、选择优良个体、优良个体配
% 对、通过交叉变异产生下一代个体，不断迭代直到满足进化终止条件。
% 下方通过程序说明遗传算法进行过程，程序可见文件：Genetic_algorithm.m
clear,clc
% 参数初始化
f=@(x)1/12*sum(x.^4)-1/2*sum(x.^2)+sum(abs(x-4));% 目标函数
lb=-1e2*ones(1,6);ub=1e2*ones(1,6);%变量范围
n=length(lb);% 变量个数
% 遗传算法计算过程
mumb=50;% 种群大小
mumb1=10;% 优异数量
mn=1e3;% 最大迭代次数
cross0=3;% 交叉基因个数
mumb2=20;% 交叉数
mumb3=20;% 变异数
gen0=3;% 变异基因个数
gama0=1;% 变异程度迭代变化参数
x0=zeros(mumb,n);% 种群初始化
mn1=10;% 遗传算法计算次数
result1=zeros(mumb,1);% 计算适应度
result_final=[];% 最终结果
for i=1:mn1
    tic
    for i1=1:mumb
        x0(i1,:)=lb+(ub-lb).*rand(1,n);
    end
    x1=x0;% x1用来迭代计算
    for i1=1:mumb
        result1(i1)=f(x1(i1,:));
    end
    gama=gama0;
    for i1=1:mn % 迭代计算
        [result3,k]=sort(result1);
        xk=x1(k,:);
        xk1=xk(1:mumb1,:);% 找到优异个体
        % 交叉操作
        xk0=zeros(mumb2+mumb3,n);
        for i2=1:mumb2
            m1=floor(mumb1*rand()+1);
            m2=floor(mumb1*rand()+1);
            xk2=xk1(m1,:);xk3=xk1(m2,:);
            for i3=1:cross0 % 交叉过程
                m3=floor(n*rand()+1);
                temp=xk2(m3);
                xk2(m3)=xk3(m3);
                xk3(m3)=temp;
            end
            xk0(i2,:)=xk2;
        end
        % 变异操作，此处可保持变异程度不变，或者令变异程度随迭代增加减小
        for i2=1:mumb3
            m1=floor(mumb1*rand()+1);
            xk0(mumb2+i2,:)=xk1(m1,:);
            for i3=1:gen0
                m2=floor(n*rand()+1);
                z=-1+2*rand();
                temp=xk1(m1,m2)+z*gama*(ub(m2)-lb(m2))/4;
                if(temp>ub(m2))
                    temp=xk1(m1,m2)+1/2*(ub(m2)-xk1(m1,m2));
                elseif(temp<lb(m2))
                    temp=xk1(m1,m2)+1/2*(lb(m2)-xk1(m1,m2));
                end
                xk0(mumb2+i2,m2)=temp;
            end
        end
        x1=[xk1;xk0];
        for i2=1:mumb
            result1(i2)=f(x1(i2,:));
        end
        gama=gama*gama;
    end
    toc
    result_final=[result_final;[result1,x1]];
end
% 可见遗传算法每次计算有较大概率准确的找到全局较优解所在位置，由于较优个体需
% 要变异操作产生，则变异程度尽量不随迭代次数变化。
% 以上程序中交叉和变异操作在原始变量基础上进行，在解决其他问题时，也可以将变
% 量转化为二进制进行交叉变异，不过计算量将会更大。
%% 3.粒子群算法Particle swarm optimization
% 粒子群算法是一种进化计算方法，源于对鸟群捕食行为研究。粒子群算法一般包括以
% 下参数：种群个数，可行解维数，最大迭代次数，位置参数限制，速度限制，惯性权
% 重，自我学习因子，群体学习因子。
% 下方通过程序说明粒子群算法进行过程，程序可见文件：Particle_swarm_optimization.m
clear,clc
% 参数初始化
f=@(x)1/12*sum(x.^4)-1/2*sum(x.^2)+sum(abs(x-4));% 目标函数
lb=-1e2*ones(1,6);ub=1e2*ones(1,6);%变量范围
n=length(lb);% 变量个数
% 粒子群算法计算过程
N=100;% 粒子个数
mn=10;% 求解次数
mn1=1e3;% 最大迭代次数
Vm=10;% 最大速度
c1=2;% 自我学习因子
c2=2;% 群体学习因子
w=1;% 惯性权重
x0=zeros(N,n);x1=x0;% 粒子群初始化
v=zeros(N,n);% 速度初始化
result0=zeros(N,1);% 储存迭代结果
result_final=[];% 储存最后结果
for i=1:mn
    tic
    for i1=1:N
        x1(i1,:)=lb+rand(1,n).*(ub-lb);
    end
    pbest=x1;% 粒子个体最佳位置初始化
    for i1=1:N
        result0(i1)=f(x1(i1,:));
    end
    [result1,k]=sort(result0);
    gbest=x1(k(1),:);% 群体最佳位置初始化
    gb0=f(gbest);
    for i1=1:mn1
        % 速度与位置迭代
        for i2=1:N
            v(i2,:)=w*v(i2,:)+c1*rand()*(pbest(i2,:)-x1(i2,:))+...
                c2*rand()*(gbest-x1(i2,:));
            if(norm(v(i2,:))>Vm)
                v(i2,:)=v(i2,:)/norm(v(i2,:))*Vm;
            end
        end
        x1=x1+v;
        % 更新个体最优位置
        for i2=1:N
            pb0=f(x1(i2,:));
            if(pb0<result0(i2))
                pbest(i2,:)=x1(i2,:);
                result0(i2)=pb0;
            end
        end
        % 更新群体最优位置
        [result1,k]=sort(result0);
        gbest1=x1(k(1),:);
        if(result1(1)<gb0)
            gbest=gbest1;
            gb0=result1(1);
        end
    end
    result_final=[result_final;[gb0,gbest]];
    toc
end
% 可以看到粒子群算法可以准确，稳定的找到较优解位置。
%% 4.禁忌搜索算法Tabu search algorithm
% 禁忌搜索算法模拟人类短期记忆，算法中使用了如下的概念：禁忌表，用来存放禁忌
% 对象的表，是算法得以进行的基本前提；禁忌对象，禁忌表中被禁的自变量位置；禁
% 忌长度，禁忌表的长度。
% 下方通过程序说明禁忌搜索算法进行过程，程序可见文件：Tabu_search_algorithm.m
clear,clc
% 参数初始化
f=@(x)1/12*sum(x.^4)-1/2*sum(x.^2)+sum(abs(x-4));% 目标函数
lb=-1e2*ones(1,6);ub=1e2*ones(1,6);%变量范围
n=length(lb);% 变量个数
% 禁忌搜索算法计算过程
mn=10;% 计算次数
mn1=1e3;% 迭代次数
t0=6;% 禁忌表长度
r=(ub-lb)/10;% 对于连续型问题，设置领域范围
r1=mean((ub-lb)/1000);% 禁忌范围
tabu=zeros(t0,n);% 禁忌表
k=50;x2=zeros(k,n);f2=zeros(k,1);
result_final=[];
for i=1:mn
    tic
    x1=lb+rand(1,n).*(ub-lb);% 初始解
    f0=f(x1);
    s0=f0;% 最优值
    xbest=x1;% 最优解
    for i1=1:mn1
        for i2=1:k
            t1=0;
            while t1==0
                z=-1+2*rand(1,n);
                xk=x1+z.*r;
                t2=0;
                for i3=1:t0
                    if(sum((xk-tabu(i3,:)).^2)>r1)
                        t2=t2+1;
                    end
                end
                if(t2==t0)
                    t1=1;
                end
            end
            x2(i2,:)=xk;
        end
        for i2=1:k
            f2(i2)=f(x2(i2,:));
        end
        [f3,k1]=sort(f2);
        tabu1=tabu(1:end-1,:);
        tabu=[x2(k1(1),:);tabu1];
        if(f3(1)<s0)
            s0=f3(1);
            xbest=x2(k1(1),:);
        end
        x1=x2(k1(1),:);
        f0=f3(1);
    end
    result_final=[result_final;[s0,xbest]];
    toc
end
% 可以看到禁忌搜索算法能够找到全局最优解所在区域。
%% 5.蚁群算法Ant colony algorithm
% 蚁群算法是一种模拟蚂蚁觅食行为的优化算法，主要用于解决旅行商(TSP)类问题。
% 算法基本过程为：在每个城市随机放置蚂蚁，蚂蚁在路径释放信息素，碰到还未走过
% 的路口，随机选择路线，并释放与路径长度成反比的信息素。
% 本节代码参考选用CSDN博主「LXzerorb」的原创文章
% 原文链接如下：https://blog.csdn.net/liuxin0108/article/details/115936358
% TSP问题代码可见文件：Ant_colony_algorithm.m
clear,clc
% 初始化参数
m=50;% 蚂蚁个数
Alpha=3;% 信息素重要程度参数
Beta=3;% 启发式因子重要程度参数
Rho=0.1;% 信息素蒸发系数
G=200;% 最大迭代次数
Q=100;% 信息素增加强度系数
load city.mat
C=city;% 城市坐标
n=size(C,1);% n表示问题的规模(城市个数)
D=zeros(n);% D表示两个城市距离间隔矩阵
for i = 1:n
    for j = 1:n
        if i ~= j
            D(i,j) = ((C(i,1)-C(j,1))^2+(C(i,2)-C(j,2))^2)^0.5;
        else
            D(i,j) = eps;
        end
        D(j,i) = D(i,j);
    end
end
Eta=1./D;% Eta为启发因子，这里设为距离的倒数
Tau=ones(n);% Tau为信息素矩阵
Tabu=zeros(m,n);% 存储并记录路径的生成
NC=1; %迭代计数器
R_best=zeros(G,n); %各代最佳路线
L_best=inf.*ones(G,1); %各代最佳路线的长度
figure(1);% 优化解，绘制图像
% 蚁群算法迭代计算过程
while NC<=G
    % 将m只蚂蚁放到n个城市上
    Randpos=[];
    for i = 1:(ceil(m/n))
        Randpos=[Randpos,randperm(n)];
    end
    Tabu(:,1) = (Randpos(1,1:m))';
    % m只蚂蚁按概率函数选择下一座城市，完成各自的周游
    for j=2:n
        for i=1:m
            visited=Tabu(i,1:(j-1));% 已访问的城市
            J=zeros(1,(n-j+1));% 待访问的城市
            P=J;% 待访问城市的选择概率分布
            Jc=1;
            for k=1:n
                if isempty(find(visited==k, 1))
                    J(Jc)=k;
                    Jc=Jc+1;
                end
            end
            % 计算待选城市的概率分布
            for k=1:length(J)
                P(k)=(Tau(visited(end),J(k))^Alpha)...
                    *(Eta(visited(end),J(k))^Beta);
            end
            P=P/(sum(P));
            % 按概率原则选取下一个城市
            Pcum=cumsum(P);
            Select=find(Pcum >= rand);
            to_visit=J(Select(1));
            Tabu(i,j)=to_visit;
        end
    end
    if NC>=2
        Tabu(1,:)=R_best(NC-1,:);
    end
    % 记录本次迭代最佳路线
    L=zeros(m,1);
    for i=1:m
        R=Tabu(i,:);
        for j=1:(n-1)
            L(i)=L(i)+D(R(j),R(j+1));
        end
        L(i)=L(i)+D(R(1),R(n));
    end
    L_best(NC)=min(L);
    pos=find(L==L_best(NC));
    R_best(NC,:)=Tabu(pos(1),:);
    % 更新信息素
    Delta_Tau=zeros(n,n);
    for i=1:m
        for j=1:(n-1)
            Delta_Tau(Tabu(i,j),Tabu(i,j+1))= ...
                Delta_Tau(Tabu(i,j),Tabu(i,j+1))+Q/L(i);
        end
        Delta_Tau(Tabu(i,n),Tabu(i,1))= ...
            Delta_Tau(Tabu(i,n),Tabu(i,1))+Q/L(i);
    end
    Tau = (1-Rho).*Tau+Delta_Tau;
    % 禁忌表清零
    Tabu=zeros(m,n);
    % 历代最优路线
    for i=1:n-1
        plot([C(R_best(NC,i),1), C(R_best(NC,i+1),1)],...
            [C(R_best(NC,i),2), C(R_best(NC,i+1),2)],'bo-');
        hold on;
    end
    plot([C(R_best(NC,n),1), C(R_best(NC,1),1)],...
        [C(R_best(NC,n),2), C(R_best(NC,1),2)],'ro-');
    title(['优化最短距离:',num2str(L_best(NC))]);
    hold off;
    pause(0.005);
    NC = NC+1;
end
% 输出结果
Pos = find(L_best==min(L_best));
Shortest_Route=R_best(Pos(1),:);% 最佳路线
Shortest_Length=L_best(Pos(1));% 最佳路线长度
figure(2),
plot(L_best)
xlabel('迭代次数')
ylabel('目标函数值')
title('适应度进化曲线')
% 可见，相比较模拟退火算法，蚁群算法将会更快更准确找到较优解。
%% 6.简单问题优化matlab优化方法
% 本节介绍matlab内置的求解各类优化问题的函数。
% (1)求解一维函数最小值：fminbnd
% fminbnd 是一个函数文件。算法基于黄金分割搜索和抛物线插值方法。除非左端点x1
% 非常靠近右端点 x2，否则 fminbnd 从不计算 fun 在端点处的值，因此只需要为x在
% 区间x1<x<x2中定义fun。
% 如果最小值实际上出现在x1或x2处，则fminbnd返回区间(x1,x2)内部靠近极小值的点x。
clear,clc
f=@(x)x*sin(x);% 目标函数
[x,fval]=fminbnd(f,0,1e2);% 在区间(0,1e2)中寻找最小值
% (2)求解多维无约束函数最小值：fminunc,fminsearch
% x=fminunc(fun,x0)在点x0处开始并尝试求fun中描述的函数的局部最小值x。点x0可
% 以是标量、向量或矩阵。
clear,clc
f=@(x)1/12*sum(x.^4)-1/2*sum(x.^2)+sum(abs(x-4));% 目标函数
x0=4*ones(1,6);% x0初始点
options=optimoptions('fminunc','Algorithm','quasi-newton');% 此处可以选择求
% 解算法，选项是 'quasi-newton'（默认值）或 'trust-region'。默认情况下，
% quasi-newton 算法使用具有三次线搜索过程的 BFGS 拟牛顿法，'trust-region'为
% 信赖域算法,使用时需要提供梯度信息。
[x,fval]=fminunc(f,x0,options);
% fminsearch仅对实数求最小值，即向量或数组x只能由实数组成，并且f(x)必须只返回
% 实数。当x具有复数值时，将x拆分为实部和虚部。
% 使用fminsearch求解不可微分的问题或者具有不连续性的问题，尤其是在解附近没有
% 出现不连续性的情况下。fminsearch 使用 Lagarias 等的单纯形搜索法.
[x,fval]=fminsearch(f,x0);
